/*
    Copyright 2006-2012 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file
    Implementation of the functions generating the adaptive grid.
    \ingroup sfrhist */

#include "config.h"
#include <time.h>
#include <unistd.h> 
#include <sstream>
#include <fstream>
#include <string>
#include "preferences.h"
#include "sphgrid.h"
#include "grid-fits.h"
#include "CCfits/CCfits"
#include "misc.h"
#include "hpm.h"
#include "fits-utilities.h"

using namespace std;
using namespace CCfits;
using namespace mcrx;


/** This function does the makegrid initialization and creates the
    grid from the supplied snapshot. The grid creation parameters are
    set as preferences keywords and are commented in the function
    body. */
auto_ptr<nbody_data_grid> init_makegrid(Preferences& p, Snapshot& snap)
{
  // Create the grid
  const vec3d min (p.getValue("grid_min", vec3d ()));
  const vec3d max (p.getValue("grid_max", vec3d ()));
  auto_ptr<nbody_data_grid> g(new nbody_data_grid (min, max));

  // refinement criteria
  const string output_file =
    word_expand (strip_suffix (p.getValue("output_file", string ()))
		 + ".fits")[0];
  // with grid emission, the grid is refined using luminosity and the
  // SED is saved in the grid cells.  Without, the particle list is
  // kept and used for emission directly.
  const bool use_grid_emission = p.defined("use_grid_emission")?
    p.getValue("use_grid_emission", bool ()): false;
  // These are the "relative" tolerances, expressed in units of Sigma
  // of the sub sells quantities/average quantity.
  const T_float gas_tolerance (p.getValue("gas_tolerance", T_float ()));
  const T_float metal_tolerance (p.getValue("metal_tolerance", T_float ()));
  const T_float L_bol_tolerance = use_grid_emission?
    p.getValue("L_bol_tolerance", T_float ()): blitz::huge (T_float ());
  const nbody_data_cell tolerance (gas_tolerance, vec3d(0,0,0), L_bol_tolerance, 
				   0, metal_tolerance, 0, vec3d(0,0,0), 0, 0, 0, 0, 0);
  // quantities needed for the absolutes tolerance calculation.  we
  // read the number of rays as a floating point so we can make it
  // ridiculously large.  If any of the keywords are not defined, it
  // has a default value meaning that the test is disabled and is NOT
  // used to unify any cells
  const T_float n_rays = p.defined("n_rays_estimated")?
    p.getValue("n_rays_estimated", T_float ()): blitz::huge (T_float ());
  // This is the opacity of *metals* now.  The opacity of gas is
  // controlled by this and the gas_metallicity.  Essentially, the
  // opacity for gas will be the opacity for metals multiplied by the
  // gas metallicity. Note that setting the metallicity here only
  // affects the grid refinement, it does not imply that gas will have
  // opacity during the radiated transfer calculation which is set in
  // the mcrx configuration file.
  const T_float kappa = p.defined("opacity" )?
    p.getValue("opacity", T_float ()): blitz::huge (T_float ());
  const T_float gas_metallicity = p.defined("gas_metallicity" )?
    p.getValue("gas_metallicity", T_float ()): 0;
  // Calculate V_grid ^ 2/3
  const T_float vg23 = pow (product (vec3d (g->getmax() - g->getmin())), 2./3);
  const nbody_data_cell absolute
    (vg23/(n_rays*kappa*gas_metallicity), vec3d(0,0,0),
     use_grid_emission? 1./n_rays: blitz::huge (T_float ()),
     0, vg23/(n_rays*kappa), 0, vec3d(0,0,0), 0, 0, 0, 0, 0);
  // This specifies what optical depth the grid should resolve.
  const T_float tau_tolerance = p.defined("tau_tolerance") ?
    p.getValue("tau_tolerance", T_float()) : blitz::huge(T_float());
  const nbody_data_cell tol_quantity(tau_tolerance/(kappa*gas_metallicity),
				     vec3d(0,0,0),
				     blitz::huge(T_float()),0,
				     tau_tolerance/kappa,
				     0,vec3d(0,0,0), 0,0,0,0,0);
  const T_float size_factor (p.getValue("size_factor", T_float ()));
  const int max_level (p.getValue("max_level", int ()));
  if (p.defined("n_threads"))
    g->use_threads(p.getValue("n_threads", int ()));
  if (p.defined("work_chunk_levels"))
    g->set_work_chunk(p.getValue("work_chunk_levels", int ()));
  if (p.defined("mcrx_data_dir"))
    particle_base::inp_dir = 
      word_expand (p.getValue("mcrx_data_dir", string ())) [0];

  // open output file and output the makegrid-specific info
  // Note that we now open with ! so any previous file is overwritten
  auto_ptr<FITS> output (new FITS("!"+output_file, Write));
  ExtHDU* makegrid_HDU = output->addTable("MAKEGRID", 0);
  makegrid_HDU->writeDate();

  // tolerance is scale-free, it's the fractional deviation allowed.
  // Abs_tol is not.  The L_bol member is essentially 1/n_rays, and it
  // should be multiplied by the total L_bol in which case it will
  // mean the luminosity deviation that will cause one ray to be
  // emitted.  The m_g member is trickier.  In essence, we want to
  // ignore changes that are so small that they will not affect one
  // interaction, but that depends not only on the local opacity but
  // also on the radiation intensity in the cell, something we clearly
  // do not know. If we approximate the radiation field to be uniform
  // and isotropic (clearly wrong), the number of rays traversing a
  // cell i should be n_i = n_rays*(V_i/V)^2/3.  An optical depth that
  // means that one of these will scatter will be tau = 1/n_i, so a
  // smaller deviation in this is negligible.  The problem is now that
  // this depends on the dust opacity and dust to gas ratio, so we
  // will have to supply these.  I think the end result is that if the
  // mass deviation is <V^2/3/(n_rays*kappa), where V is the total
  // grid volume and kappa is opacity per unit density.  Interestingly
  // enough, the result ends up being independent of the volume of the
  // cell itself.

  // make a tolerance checker object
  tolerance_checker tc(tolerance, absolute, gas_metallicity,
		       tau_tolerance/kappa);

  // now go ahead and load the snapshot. This used to load from the
  // file but now loads from the snapshot object
  g->load_snapshot(snap, max_level, size_factor, tc,
		   makegrid_HDU, use_grid_emission);

  // and now save
  cout << "Saving grid to file " << output_file << endl;
  g->save_structure(*output, "GRIDSTRUCTURE", p.getValue("save_cellcodes", false, true));
  // save length units in grid structure HDU
  g->save_data(*output, "GRIDDATA", "INTEGRATED_QUANTITIES");
  output.reset(0);

  // optionally, dump the grid structure
  if (p.defined("dump_grid_structure")) {
    ofstream o (word_expand (p.getValue("dump_grid_structure", 
					string ()))[0].c_str());
    g->structure_plot(o);
  }

  cout << "\nMakegrid complete.\n" << endl;
  return g;
}


  
